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ABSTRACT 


Changes  in  the  globally  integrated  absolute  angular  momentum  of  the  atmosphere 
were  computed  from  the  Fleet  Numerical  Oceanography  Center  NOGAPS  wind  analyses 
and  compared  to  astronomically  measured  changes  in  length  of  day  ( LOD )  obtained 
from  the  L'.S.  Naval  Observatory,  Washington  D.C.  The  two  time  series  were  subjected 
to  both  time  and  frequency  domain  analysis.  In  the  time  domain,  digital  filters  were  used 
to  isolate  seasonal  and  subseasonal  components.  In  the  frequency  domain,  energy  den¬ 
sity,  coherence  and  phase  were  computed  over  periods  from  2  days  to  1000  days.  Over 
90%  of  the  total  variance  in  astronomically  determined  LOD  can  be  explained  by 
meteorological  phenomena.  Fluctuations  in  LOD  are  coherent  and  in  phase  with  fluc¬ 
tuations  in  the  globally  integrated  angular  momentum  of  the  Earth's  shell  (crust,  mantle 
and  oceans;  liquid  core  is  excluded)  at  almost  all  periods  less  than  365  days.  Annual 
fluctuations  in  LOD  appear  to  originate  in  the  midlatitudes  and  propagate  equatorward. 
Subseasonal  fluctuations  (30  to  100  day  periods)  appear  to  be  a  tropical  phenomena. 
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I.  INTRODUCTION  AND  THEORY 


A.  BACKGROUND 

Fluctuations  in  the  Earth's  length  of  day  (LOD)  were  surmised  by  early  astronomers 
who  noted  cumulative  irregularities  in  the  Earth's  expected  latitudinal  position  relative 
to  daily  star  transits.  Newtonian  physics  explained  only  the  tidal  component  of  these 
short  period  changes  in  the  Earth's  angular  momentum.  In  the  late  1930's,  pendulum 
clocks  were  accurate  enough  to  measure  daily  star  transits  to  within  1  ms.i  Stoyko  (1937, 
cited  in  Munk  and  MacDonald,  1960)  reported  a  2  ms  cumulative  seasonal  variation  in 
LOD.  He  noted  that  the  length  of  day  measured  in  July  was  longer  than  that  measured 
in  May  by  2  ms.  Atomic  clocks,  developed  in  the  1950's,  improved  time  keeping  accu¬ 
racy  by  several  orders  of  magnitude.  Present  day  cesium  clocks  have  accuracies  better 
than  10-*  s  per  day  (or  .00001  ms)  (Ramsey,  1988).  At  the  same  time,  improvements  in 
astronomical  measurement  methods  led  to  more  precise  determination  of  star  and  quasar 
transits. 

A  number  of  studies  have  been  undertaken  to  link  the  observed  changes  in  LOD  to 
meteorological  phenomena.  Starr  (1948,  cited  in  Munk  and  MacDonald,  1960)  pro¬ 
posed  that  the  observed  changes  in  the  Earth's  LOD  occurred  because  the  solid  Earth 
and  atmosphere  exchanged  angular  momentum  on  a  seasonal  basis.  Madden  and  Julian 
(1971)  found  40-50  day  oscillations  in  pressure  and  zonal  winds  centered  in  the  tropics, 
and  Lambeck  and  Cazenave  (1974)  subsequently  linked  these  atmospheric  oscillations 
to  changes  in  LOD.  Lambeck  (1980)  identified  the  amplitudes  and  phases  of  the  seasonal 
changes  and  listed  other,  non  meteorological  components  which  could  be  responsible  for 
LOD  changes.  Rosen  and  Salstein  (1983)  correlated  short  period  (less  than  seasonal) 
hemispherical  changes  of  atmospheric  angular  momentum  and  LOD  with  polar  and 
subtropical  jet  streams.  Eubanks  et  al.  (1985)  performed  a  spectral  analysis  of  atmo¬ 
spherically  derived  angular  momentum  with  a  composite  of  astronomic,  satellite  and 
Lunar  Laser  Ranging  (LLR)2  LOD  data.  Barnes  et  al.  (1983)  extensively  studied 

1  The  best  gravity  pendulum  clocks  have  accuracies  within  1  ms  per  day  or  one  part  in  10s 
(Ramsey,  1988).  This  is  barely  accurate  enough  to  detect  the  fairly  large  seasonal  changes  in  the 
LOD. 

2  The  astronauts  left  a  comer  reflector  on  the  moon  which  reflects  a  LASER  signal  back  for 
extremely  accurate  short  term  measurement  of  Earth  rotation  parameters. 


changes  in  LOD  over  the  period  1979-1981.  Morgan  et  al.  (1985)  used  least  squares 
analysis  to  study  the  seasonal  as  well  as  subseasonal  nature  of  oscillations  in  the  LOD 
from  1981  to  1984.  Swinbank  (1985)  examined  the  torques  which  produce  changes  in 
LOD.  Wolf  and  Smith  (1987),  and  Eubanks  et  al.  (1986)  compared  changes  in  atmo¬ 
spheric  angular  momentum  with  astronomically  determined  LOD  during  the  El  Nino  of 
1983. 

Most  studies  have  examined  atmospheric  angular  momentum  and  LOD  using  time 
measurements  from  the  Bureau  International  de  l'lleure  (Bill),  augmented  by  LLR  and 
other  precise  rotation  measurement  methods,  to  compare  with  gridded  zonal  wind  data 
from  the  U.S.  National  Meteorological  Center  (NMC)  or  the  European  Centre  for  Me¬ 
dium  Range  Weather  Forecasts  (ECMWF).  This  paper  analyzes  the  wind-derived  at¬ 
mospheric  angular  momentum  budget  for  the  period  January  1983  to  July  1987  and  the 
simultaneous  changes  in  tne  Earth's  LOD  using  time  data  from  the  U.S.  Naval  Observ¬ 
atory  (USNO)  and  zonal  wind  data  from  the  U.S.  Navy  Fleet  Numerical  Oceanography 
Center  NOGAPS  model.3  The  atmospheric  angular  momentum  changes  are  then  exam¬ 
ined  over  two  latitude  bands  to  investigate  the  horizontal  structure  of  the  seasonal  and 
subseasonal  components  of  the  atmospheric  angular  momentum  and  contributions  to 
the  LOD  . 

B.  ASTRONOMIC  LENGTH  OF  DAY 

Astronomic  length  of  day  is  measured  by  comparing  elapsed  time  for  meridional 
transits  of  stars  or  quasars  with  a  reference  time  determined  by  atomic  clocks.  The  time 
it  takes  for  exactly  one  rotation  of  the  Earth  (360  degrees  of  longitude  change)  with  re¬ 
spect  to  the  fixed,  inertial  reference  system  of  distant  stars  is  a  sidereal  day  with  a  mean 
value  of  86,164.0905  seconds.  Earth  time  is  generally  kept  with  respect  to  rotation  about 
the  sun:  the  solar  day.  A  mean  solar  day  is  slightly  longer  and  subject  to  more  variability 
than  the  sidereal  day  because  the  Earth  is  m  motion  about  the  sun.  Time  for  a  solar  day 
( LOD0 )  is  defined  as  exactly  86,400  seconds.  LOD  is  defined  as  the  actual  elapsed  time 
(measured  by  cesium  clocks)  for  a  meridional  transit  of  a  point  on  Earth  with  respect  to 
an  inertial  frame.  Since  this  is  with  respect  to  sidereal  time,  a  conversion  is  applied  to 
make  it  consistent  with  solar  time  keeping. 


3  NOGAPS  is  Naval  Operational  Global  Atmospheric  Prediction  System.  It  is  a  mathematical 
model  which  generates  a  global  wind  and  pressure  analysis  on  a  2  Vi  by  2  Vi  degree  grid  spacing. 


LOD  is  related  to  the  Earth  rotation  rate  (o>)  by: 


=  £2  • 


d(lT\) 
d(IAl)  ’ 


where  LOD  and  UT1  are  both  measures  of  the  actual  elapsed  time  of  the  solar  day,  £2 
is  the  mean  sidereal  rotation  rate  (7.2921  x  10~5  $-')  of  the  Earth,  and  I  AT  is  the  precise 
time  reference  determined  from  atomic  clocks.  671  is  an  angular  measurement  of 
Earth's  polar  rotation  with  respect  to  sidereal  time  passage  which  is  usually  converted 
to  time. 

Changes  in  the  Earth's  rotation  rate  with  respect  to  accepted  standard  values  are:-* 


<5<y  =  £2 


L0Do 


LODa 

_  o - — 

7  r\ 


d(UT\  -  7.47;, d) 


Taking  the  difference  between  LOD  and  the  standard,  LODita,  to  be  6L0D  produces: 

,  „  SLOD  „  SLOD  n  d(UTl  -  IAT,ld) 

5“~~q-wd~  a~L0i£~a - 7, - •  (3) 

where  LODM ,  is  taken  equal  to  LODa  and  I A  Tni  is  taken  equal  to  a  sidereal  day,  although 
other  conventions  could  be  used  with  the  relationship  still  holding  (Eubanks  et  al.,  1985). 
The  approximation  in  (3)  is  valid  because,  using  a  time  step  of  one  sidereal  day,  LO D0, 
the  ratio  of  SLOD  to  LOD  is  the  same  as  the  ratio  of  <5 LOD  to  LOD0  to  within  10~8s. 
Using  finite  differences  to  approximate  the  UT 1  derivative: 


SLOD  =  —  LOD, 


{UT\  -  IATsld\  -  {UTl  -  !ATnd\_Al 
At 


Eq.  4  relates  daily  observations  of  star  and  quasar  transits  directly  with  the  precise 
elapsed  time  determined  by  atomic  clocks.  It  is  computed  on  a  daily  basis  (At  =  1  day) 
from  astronomic  data  gathered  all  over  the  world  to  help  eliminate  single  site  bias.  This 

4  By  subtracting  standard  values,  which  are  the  international  definitions  of  the  solar  and 
sidereal  days,  the  remaining  quantity  represents  a  change  from  the  standard  as  well  as  a  change  with 
respect  to  exactly  one  solar  day.  The  actual  measured  length  of  day  is  taken  with  reference  to  a 
beginning  epoch  (l  January  1900,  by  convention).  Thus  the  actual  length  of  day  (LOD)  is  the 
running  sum  of  the  changes  in  LOD  from  the  beginning  of  the  epoch.  Thus,  the  term  SLOD,  as 
developed  in  Eq.  4  represents  a  time  derivative  only  in  the  sense  of  it  being  the  first  difference  of 
the  LOD  time  series.  If  the  first  difference  were  not  used,  the  numbers  associated  with  the  length 
of  day  would  be  so  large  as  to  be  cumbersome. 
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quantity  is  compared  with  the  atmospherically  derived  equivalent,  SLOD,lm,  developed 
in  the  next  section. 

C.  ATMOSPHERIC  ANGULAR  MOMENTUM  AND  LENGTH  OF  DAY 

Over  long  periods  of  time  (centuries  and  longer)  Earth  angular  momentum  (.UE„lh) 
slowly  transfers  to  the  moon  through  solid  body  tidal  friction;  hence  LOD  increases 
slowly  and  ,\/Ejrlh  decreases.  For  periods  less  than  10  years,  this  effect  is  negligible;  .UE,nh 
is  essentially  conserved  (Eubanks  et  al.\  1985,  Lambeck,  1980).  Thus,  changes  in  ,V/Earth, 
denoted  by  (5d/Eartn  will  sum  to  zero: 

<5  -V4arth  =  ^shell  +  ^ocean  +  +  &  ^core  =  0.  (5) 

Since  5Mmn  and  S interact  in  an  oscillatory  manner  with  a  period  of  10  years  or 
greater,  <5  .!/«,„  appears  in  short  data  records  (under  five  years)  as  a  linear  trend  and  is 
easily  removed  by  detrending.  The  three  remaining  momentum  elements:  c5,\/llm,  5  A/lfiell 
and  d.V/oeMn  interact  to  produce  changes  in  the  Earth  measured  LOD.  is  difficult 

to  measure  directly.  To  exchange  momentum  with  the  rotating  Earth,  the  oceans  must 
have  unimpeded  zonal  flow  about  the  globe.  Only  the  Antarctic  Circumpolar  Current 
(ACC)  has  such  characteristics.  Lambeck  (1980)  agrees  with  the  earlier  findings  of  Munk 
and  MacDonald  (1960)  that,  at  most,  the  ACC  could  have  a  semiannual  amplitude  5% 
of  the  atmospheric  amplitudes.  Thus,  with  <5A/E,„h  =  0,  and  disregarding  inputs  from 
the  ocean  and  core: 

<^shell  =  -  •  (6) 


The  atmosphere  of  the  Earth  represents  about  ICE6  of  the  total  mass  of  the  solid 
Earth  (solid  Earth  includes  crust,  mantle,  oceans  and  atmosphere;  everything  except  the 
liquid  core).  Because  of  the  large  difference  in  atmospheric  and  Earth  masses,  relatively 
large  changes  in  the  zonal  velocity  field  of  the  atmosphere  are  required  to  produce  small 
changes  in  M,mu. 

The  angular  velocity  (co)  of  the  Earth's  shell  (solid  Earth  minus  the  atmosphere)  is: 


M, 


co  =  ■ 


shell 


'shell 


(7) 


where  is  the  angular  momentum  of  the  Earth's  shell  and  7lhell  is  the  polar  moment 
of  inertia.  Again,  after  subtracting  standard  values  from  <5A/,h„|,  changes  in  co  are: 


(8) 


Ou)  = 


Overbars  indicate  average  values  with  respect  to  time.  Since  7(hell  changes  little  on  time 
scales  less  than  decades,  it  is  taken  as  a  constant:  7.103  x  I0r  kg  Substituting 
-  from  (6)  into  (8)  and  combining  with  (3)  yields: 


SLOD  = 


L00o  S\tMm 


The  assumptions  leading  to  Eq.  (9)  imply  that  fluctuations  in  the  length  of  day  are  di¬ 
rectly  proportional  to  fluctuations  in  the  atmospheric  angular  momentum.  One  purpose 
of  this  thesis  is  to  see  if  independent  observations  of  SLOD  and  are  indeed  related 
as  predicted  by  (9). 

The  atmospheric  angular  momentum  at  any  time  can  be  measured  by  evaluating  a 
volume  integral  of  the  east-west  (u)  component  of  the  wind  from  the  surface  of  the  Earth 
to  the  top  of  the  atmosphere. 


Matm  =  |02ffJ_\|rr_7pU(r’  4>'  d)  ^  QOs2{4>)  dr  d<t>  d6'  (10) 

where  r  is  distance  from  the  center  of  the  Earth  and  U  is  the  absolute  zonal  velocity, 
U  =  rQ.  cos  <t>  +  u,  and  u  the  zonal  wind.  The  hydrostatic  relationship  allows  pressure 
coordinates  to  be  substituted  for  p  dr: 


dr  =  .  (1.) 

Since  the  depth  of  the  atmosphere  is  very  small  compared  to  the  radius  of  the  Earth,  r 
is  accurately  approximated  by  the  mean  radius  of  the  Earth,  RE  (6.37  x  10*  m).  With  the 
above  substitutions  into  (10),  :V/tun  becomes: 


Bl 


Matm  =  J  J02*J_  COS(<*>)  +  »>Rl  COS2 (4))  dp  d*  dO  .  (12) 


Separating  the  two  parts  of  the  integral  and  simplifying, 


2nRl  fiooombfY  ■> 

A4lm  =  -^j100mb  cotmuWdp 

4  n 

+  Q  —jr  cos\<t>)  A  P  d<j>  dd.  ( 1 3) 

The  first  term  is  the  wind  or  motion  term  (Barnes  et  al.  1983).  The  zonal  wind  av¬ 
eraged  by  latitude  bands  is  indicated  by  [u]  .  Integration  limits  of  100  mb  and  1000  mb 
for  the  motion  term  were  selected  because  the  NOGAPS  wind  analysis  extends  to  only 
100  mb  and  because  use  of  surface  winds  would  add  unnecessary  complexity  to  the  at¬ 
mospheric  angular  momentum  calculations.  Since  winds  above  100  mb  constitute  less 
than  10%  of  the  mass  of  the  atmosphere,  exclusion  of  the  stratospheric  winds  should 
introduce  fairly  small  errors.  Rosen  and  Salstein  (1985)  found  the  stratospheric  contrib¬ 
ution  to  to  be  positive  at  semiannual  periods  and  negative  at  annual  periods.  In 
other  words,  the  stratosphere  opposes  the  troposphere  (180°  out  of  phase)  at  annual 
cycles,  and  adds  to  tropospheric  angular  momentum  at  semiannual  cycles.  The  ampli¬ 
tudes  of  these  signals  are  less  than  10%  of  the  tropospheric  term  for  annual  periods,  and 
30%  of  the  tropospheric  term  for  the  semiannual  term.  Surface  winds  are  not  used  be¬ 
cause  of  conflicts  where  surface  pressure  is  near  1000  mb  (most  of  the  globe).  In  that 
case  the  near  surface  winds  would  be  doubly  represented. 

The  second  term  of  (13)  is  the  pressure,  or  matter,  term  (Barnes  et  al.,  1983).  Most 
studies  to  date  have  dropped  the  pressure  term  for  lack  of  accurate  surface  pressure  data. 
Barnes  et  al.  (1983)  used  the  special  observations  from  the  GARP  (Global  Atmospheric 
Research  Program)  of  1979  to  evaluate  the  pressure  term  and  its  contribution  to  Mtm. 
They  concluded  that  pressure  variations  in  Mtm  are  in  phase  with  the  wind  term  and 
fairly  constant  at  10%  of  the  wind  term  amplitude.  Indirect  studies  by  Morgan  et  al. 


(1985)  show  the  pressure  term  to  be  10%  of  the  wind  term  at  seasonal  periods  and  no 
more  than  25%  at  shorter  periods. 

Standard  values  are  subtracted  from  Afllm  to  give  <5 A/,lm  .  The  standard  value  for  both 
the  wind  and  pressure  terms  are  taken  to  be  zero.  This  is  consistent  with  the  atmosphere 
having  a  net  angular  momentum  of  zero  with  respect  to  the  Earth.  Neglecting  the 
pressure  term  in  (13)  and  using  (9),  changes  in  the  Earth  observed  angular  momentum 
caused  by  the  atmosphere  are: 


SLOD^lm 


LOD0  2 nRj  j*  lOOOmb  j.-j- 
«  /shell  g  J,00mb  J-f 


[«]  COSi(4>)  d(j)  dp. 


(14) 


Eq.  14  represents  the  change,  in  units  of  time,  that  changes  in  the  zonal  wind  will  have 
on  the  rotation  of  the  Earth.  Thus,  the  independently  determined  measures  of  5LOD 
from  (4)  and  (14)  should  be  the  same  if  the  atmosphere  is,  in  fact,  responsible  for  forcing 
short  period  changes  in  Earth  observed  length  of  day.  The  data  processing  and  the 
variance  computations  are  described  in  the  next  chapter,  and  all  of  the  results,  impli¬ 
cations  and  interpretations  are  given  in  Chapter  III. 
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II.  DATA  ACQUISITION  AND  PROCESSING 


A.  U.S.  NAVAL  OBSERVATORY  LOD  MEASUREMENTS 

Measures  of  SLOD  for  the  period  January  1983  to  July  1987  were  obtained  from  the 
U.S.  Naval  Observatory,  Washington,  D.C.  (USN'O).  Included  in  the  data  set  were  both 
the  Naval  Observatory  and  the  French  Bureau  International  de  l'Heure  (BIH)  estimates 
of  SLOD.  For  notation  purposes,  <5 LODvs>i0  will  refer  to  the  USN'O  data  and  SLODBlH 
to  the  French  data. 

B.  NOGAPS  WIND  DATA 

Twice  daily  analyses  of  the  global  zonal  winds  for  the  January  1983  through  August 
1987  observation  period  were  obtained  from  the  NOGAPS  model  at  the  U.S.  Navy  Fleet 
Numerical  Oceanography  Center,  Monterey  CA.  The  data  were  in  the  Navy  Environ¬ 
mental  Data  Network  (NEDN)  format.  Reporting  levels  for  the  winds  were  standard 
levels  plus  925  and  250mb.5  The  2  Vi  0  by  2  Zi  0  grid  produces  10,512  reporting  points 
for  the  winds.  A  FORTRAN  program  was  written  to  decode  the  NEDN  format,  ma¬ 
nipulate  the  11  by  10,512  element  matrix  and  reduce  the  data  to  a  single  value  of 
SLOD,im.  Eq.  14  was  approximated  by  zonally  averaging  the  matrix,  integrating  verti¬ 
cally  (using  a  trapezoidal  approximation  for  the  unevenly  spaced  pressure  levels),  ap¬ 
plying  the  cos2  latitude  terms  and  summing  over  latitude.  Since  the  integration  over 
latitude  was  performed  last,  the  program  was  easily  modified  to  get  SLOD,m  values  for 
specific  latitude  bands  as  well  as  the  globe  as  a  whole.  Appendix  A  contains  the  source 
code  for  the  program. 

The  archived  data  had  few  (less  than  4%)  missing  data  points.  All  of  the  missing 
observations  in  the  wind  data  were  interspersed  such  that  linear  extrapolation  of  the 
missing  winds  could  be  made.  The  resulting  time  series  was  smoothed  with  a  three-day 
moving  average  filter  and  decimated  to  1622  elements. 

Similar  steps  were  taken  in  determining  SLODtxm  values  for  latitude  bands  within 
15°  and  20°  of  the  equator,  respectively.  Values  for  the  remainder  of  the  globe  were  ob¬ 
tained  by  subtracting  the  latitude  band  series  from  the  original  6LODtm  series.  The 

5  Standard  levels  include  surface,  1000,  850,  700,  500,  400,  300,  200,  150  and  100  mb.  Surface 
winds  were  not  used  to  calculate  SLODttm  because  of  the  possible  confusion  with  the  1000  inb 
winds. 


purpose  of  the  regional  analysis  was  to  investigate  the  possible  physical  mechanisms 
contributing  to  the  LOD  and  \fam  relationships. 

C.  ANALYSIS  TECHNIQUES 

Fig.  1  on  page  22  shows  the  global  SLOD,{m  and  the  dLOD^so  plotted  together.  The 
time  data  have  an  obvious  decreasing  linear  trend,  while  the  atmospheric  data  do  not. 
The  means  of  both  series  were  removed  prior  to  plotting.  The  trend  in  the  time  data  is 
from  decade  fluctuations  in  the  LOD  caused  by  core, 'shell  coupling  (Munk  and 
MacDonald,  1960  and  Lambeck,  1980).  The  trend  was  removed  by  subjecting  both  series 
to  a  least  squares  algorithm  found  in  Bendat  and  Piersol  (1971,  equation  9.20).  The  two 
detrended  series  are  plotted  in  Fig.  2  on  page  22.  All  of  the  series  developed  in  the 
analysis  were  detrended  with  this  scheme. 

1.  Time  Domain  Analysis 

Both  series  were  subjected  to  linear  filtering  to  isolate  seasonal  and  subseasonal 
components.  First,  the  high  frequency  noise  and  tidal  oscillations  were  removed  with  a 
sine  squared  weighted  moving  average  filter  (also  known  as  a  Hanning  filter)  of  23  points 
(Fig.  3  on  page  23).  Weights  for  such  a  filter  are  determined  by  Robinson  and  Silvia, 
(1978): 

Wn  =  sin2(^y)  n  =  1,2,3 M,  (15) 

where  IV „  are  the  filter  weights  and  M  is  the  total  number  of  weights.  The  sinJ  weights 
minimize  spectral  distortion  through  unwanted  side  lobe  leakage  at  the  expense  of  a 
slightly  smeared  frequency  spectrum  (Robinson  and  Silvia,  1978).  By  keeping  the  filter 
symmetric,  phase  information  of  both  series  is  preserved.  Symmetric  application  of  the 
filter  results  in  the  general  expression  of  the  transformed  series  (Robinson  and  Silvia, 
1978,  Eq.  6.4): 

si- 1 
2 

%  =  m  _  i  ^  ^  xt-i  •  (16) 


The  frequency  response,  or  convolution,  of  a  symmetric,  moving  average  filter  of  M 
points  is  given  by: 


(17) 


H 00  =  ~~  -  (/to  +  ^/i,  cos(2«  nf)), 

i=  i 


A/  +  1 

where  the  series  h  represents  the  first  (n  =  - — - — )  weights  of  the  symmetric  series  IV„ 
and  /f„  is  the  center  weight  of  the  full  series  \t.  The  term  /ia  will  always  equal  one,  since 
it  is  the  center  weight.  Subtraction  moving  average  filters  have  a  transfer  function  of 
I  —  H(/)  and  serve  to  pass  high  frequencies.  Sine  squared  weights  produce  transfer 
functions  almost  devoid  of  side  lobes.  Plots  of  H (/)  for  each  of  the  filters  used  are  con¬ 
tained  in  Appendix  C. 

One  drawback  to  using  convolution  techniques  is  the  loss  of  A/  —  1  data  points 
(half  at  each  end  of  the  series)  with  each  pass  of  the  filter.  This  is  particularity  noticeable 
when  using  a  long  string  of  filter  weights  often  required  for  lowpass  systems.  Nonethe¬ 
less,  the  record  length  of  the  data  (1622  points),  affords  sufficient  length  to  use  fairly 
long  series  of  weighted  moving  average  filters. 

To  isolate  the  seasonal  terms  in  the  <5UODlsno  and  5LODum  series,  a  75-point 
sin2  weighted  lowpass  filter  was  applied  twice  to  both  series.  The  effect  of  successive  fil¬ 
terings  in  the  time  domain  is  multiplication  of  the  transfer  functions  in  the  frequency 
domain.  The  transition  band  with  a  double  pass  is  narrower  than  that  of  a  single  pass. 
Periods  less  than  50  days  are  completely  attenuated.  The  amplitude  of  annual  periods 
are  attenuated  by  5%  and  the  semiannual  terms  are  diminished  by  18%.  The  resulting 
series  are  plotted  in  Fig.  4  on  page  23  and  the  transfer  function  is  plotted  in  Fig.  18  on 
page  47.  To  further  separate  the  seasonal  terms  into  annual  and  semiannual  compo¬ 
nents,  both  series  were  subjected  to  a  183-point  unweighted  moving  average  filter.  The 
sin2  weights  w'ould  have  required  too  many  points  to  adequately  filter  out  the  semiannual 
terms.  Using  the  "boxcar"  unweighted  filter  of  183  points  put  the  first  zero  at  exactly  the 
semiannual  period,  allowing  almost  complete  attenuation  of  those  signals.  Disadvan¬ 
tages  of  this  type  of  filter  are  the  fairly  large  side  lobes  at  frequencies  above  the  first  zero 
crossing  (183  days).  Since  the  boxcar  has  a  fairly  narrow  transition  band,  the  remaining 
annual  terms  were  still  about  65%  of  their  original  amplitude.  The  filter  response  is 
contained  in  Fig.  19  on  page  48  in  Appendix  C,  and  the  filtered  series  is  Fig.  5  on  page 
24.  To  isolate  the  semiannual  terms,  the  series  obtained  from  the  183-point  boxcar  filter 
was  subtracted  from  the  twice  run  75-point  weighted  series.  The  resulting  semiannual 
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series  was  attenuated  less  than  15%  at  periods  of  200  days.  Response  for  this  series  is 
contained  in  Fig.  20  on  page  48. 

Oscillations  with  periods  from  20  to  100  (subseasonal)  days  were  isolated  with 
high  and  lowpass  filters.  The  highpass  filter  consisted  of  two  passes  of  a  subtracted 
moving  average  of  75  points.  Lowpass  was  achieved  with  a  single  run  of  a  23-point 
weighted  moving  average  filter.  The  resulting  series  are  plotted  in  Fig.  7  on  page  25. 
l  ig.  21  in  Appendix  C  contains  the  frequency  response  of  the  filter. 

For  latitude  band  comparisons,  the  SLODu%so  series  was  separated  into  rela¬ 
tively  high  frequency  (subseasonal)  and  relatively  low  frequency  (seasonal)  components. 
The  transition  band  was  from  1/60  day-1  to  1/150  day-'.  These  two  series  (heavy  lines 
in  Figs.  4  and  7),  5LO/)usso.„,„„„  and  5LODLSN0.,gt),M1,  were  then  individually  compared 
to  the  corresponding  6LOD,m  time  series  generated  from  tropical  and  extratropical  zonal 
winds.  These  plots  are  contained  in  Figs.  8  through  11.  Naming  conventions  for  these 
series  are: 


•  <5 LODViS0.mnnl:  Seasonal  (150  to  365  day  periods)  values  of  USNO  determined 
SLOD  . 

•  ^LODLSN0.sub„t>:  Subseasonal  (20  to  100  day  periods)  values  of  USNO  determined 
5LOD. 

•  SLOD,mls:  Atmospherically  derived  SLOD  from  15°/V  to  15°S. 

•  SLOD,tm.s:  Atmospherically  derived  SLOD  from  15°  N  to  the  North  Pole  and  15°S 
to  the  South  Pole. 

•  SLOD,m 20  and  5LOD,lm,„:  Atmospherically  derived  SLOD  with  cutofF  at 
20° ,V  and  5. 


2.  Frequency  Domain  Analysis 

Prior  to  Fourier  analysis,  each  series  was  tapered  with  a  cosine  bell  to  reduce 
leakage  at  the  ends  of  the  time  series  data.  Both  series  were  extended  to  a  length  of  2048 
points  by  the  addition  of  zeros  to  each  end  of  the  record.  The  effects  of  the  tapering  were 
removed  after  the  transform  by  dividing  each  element  by  the  ratio  of  the  area  of  the  ta¬ 
per  (the  cosine)  to  the  rectangular  boxcar  area  replaced  by  the  taper;  a  factor  of  0.875 
(Bendat  and  Piersol,  1971). 

A  record  length  of  2048  permitted  the  use  of  the  Fast  Fourier  Transform  (FFT) 
algorithm  (Bendat  and  Piersol,  1971),  significantly  cutting  down  the  computation  time 
required  for  the  transforms.  After  transforming,  the  complex  valued  cross  spectral  am- 
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plitude  coefficients  were  determined  from  the  complex  amplitude  coefficients  of  the  time 
and  atmospheric  series  (X,M  and  Y,ni).  Cross  spectral  coefficients  are: 


XY 


(n) 


2X 


(n)  T(n) 

A/ 


n  =  0,1,2..., 


(18) 


V 

where  n  represents  the  index  number  of  the  — f- 1  resolvable  frequencies  from  a  sample 
size  of  ,Y  and  X*  is  the  complex  conjugate  of  X.  The  three  complex  valued  scries, 
X,  Y,  and,  XY,  were  then  smoothed  in  order  to  increase  the  confidence  intervals  of  the 
power  spectral  density  functions  as  well  as  the  coherence  squared  (yJ)  and  phase  {<]>) 
functions.  Prior  to  smoothing,  the  spectral  estimates  were  prewhitened  by  multiplication 
by  the  square  of  the  frequency  at  each  estimate.  Jenkins  and  Watt  (1968)  state  that  es¬ 
timates  of  coherency  and  phase  should  be  based  on  white  spectra.  This  technique  (also 
used  by  Eubanks  et  ai,  1985),  compensates  for  the/'*  power  law  followed  by  both  the 
time  and  atmospheric  data.  The  complex  spectral  and  co-spectral  estimates  were 
smoothed  with  a  seven-point  unweighted  moving  average  to  give  14  degrees  of  freedom 
for  calculating  confidence  and  significance  intervals.  Energy  density  estimates  were 
formed  by  taking  the  magnitude  of  the  complex  coefficients: 
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Coherence  squared  is  defined  as: 

.2  IVV(,)I‘ 
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(20) 


Phase  is  determined  by: 
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I*n[X(n)*Y(n)] 
tan  Real[X(n)*  Y(n)]  1 


(21) 


Eqs.  19,  20  and  21  were  calculated  and  plotted  for  the  <5LO£)usno  and  3LODum  series,  as 
well  as  the  tropical  and  extratropical  3LOD,Xm  series  against  the  SLODvs>i0  .  Figs.  12 
through  16  are  the  plots  for  these  series. 
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3.  Harmonic  Analysis 

For  comparison  with  earlier  findings  of  Rosen  and  Salstein  (  1985),  Morgan  ct 
al.  (1985),  and  Langley  et  al.  (1981),  a  harmonic  analysis  was  performed  to  obtain  the 
annual  and  semiannual  terms  of  both  time  and  atmospheric  signals.  The  Fourier  anal¬ 
ysis  was  done  for  the  three  calendar  years  for  which  the  data  were  continuous:  1984, 
1985  and  1986.  Phase  and  amplitude  were  determined  from  the  sine  and  cosine  sums  of: 

■y  p365  -y  p365 

A(i)s. in  =  A-(r)  Mlnfrfdt ;  ,%cos  =  JQ  X(t)  cos(2nJ\t)dt,  (22) 

where  f  is  the  frequency  of  interest  (1/365  day1  and  1/182.5  day1)  and  t  is  the  integer 
time  in  days  from  the  beginning  to  the  end  of  the  year.  Amplitude  for  the  frequencies  is 
\A}at  +  and  phase  is  tan converted  to  days  from  1  January  of  that  year. 
Results  and  comparisons  with  earlier  findings  are  in  Table  3. 


HI.  RESULTS 


A.  GLOBAL  COMPARISON  IN  THE  TIME  DOMAIN 

The  detrended  5LODv%so  and  SLODum  series  (Figs.  2  and  3)  match  fairly  closely. 
The  linear  correlation  coefficient  for  the  two  series  with  tides  and  high  frequency  noise 
removed  is  0.990.6  From  Fig.  2  on  page  22  it  can  be  seen  that  the  signals  with  the  longer 
periods  have  larger  amplitudes.  This  is  the  time  domain  analogy  of  the  f~2  relationship 
found  in  the  spectra  of  <5LOZ\sno  and  6LOD,m.  The  lowest  period  that  correlation  can 
be  expected  is  around  20  days  because  the  23-point  filter  used  to  smooth  the  data 
completely  removes  periods  of  15  days  or  less.  Additionally,  contamination  from  the 
unwanted  tidal  signal  at  a  period  of  13  days  would  make  correlation  at  or  near  that  pe¬ 
riod  impossible. 

1.  Seasonal  Components 

Seasonal  variations  stand  out,  with  the  largest  values  of  both  SLODv^0  and 
5LODum  occurring  in  January  and  July.  Munk  and  MacDonald  (1960)  observed  these 
seasonal  elfects,  and  also  noted  that  the  amplitudes  of  the  annual  oscillations  were 
greater  in  July  than  in  January.  This  corresponds  to  the  uneven  contribution  (by  almost 
a  factor  of  two  in  Munk  and  MacDonald's  computations)  of  the  two  hemispheres  to  the 
angular  momentum  budget  during  January  and  July.  Thus,  a  typical  day  in  January  is 
almost  one  millisecond  longer  than  a  typical  day  in  July.  In  January,  strong  westerlies 
from  the  winter  (northern)  hemisphere  dominate  the  global  angular  momentum  and 
cause  the  observed  decrease  in  the  length  of  day.  In  other  words,  the  Earth  spins  slower 
because  the  atmosphere  is  spinning  faster.  In  July,  the  upper  level  westerlies  in  the 
northern  hemisphere  are  weakened  by  the  breakdown  of  the  midlatitude  baroclinic  zone. 
Meanwhile,  the  southern  hemisphere,  because  of  its  smaller  land  mass,  does  not  com¬ 
pensate  with  stronger  westerlies.  The  net  effect  is  for  a  minimum  of  atmospheric  angular 
momentum  in  July  and  a  corresponding  maximum  in  Earth  rotation  rate,  and  the  reverse 
in  January;  maximum  angular  momentum  and  minimum  earth  rotation  rate. 

6  The  linear  correlation  coefficient  is  defined  as  p  -  Cov x)lexaf.  Its  analogy  in  frequency 
space  is  coherence,  although  the  coherence  statistic  used  in  this  paper  is  more  properly  called  co¬ 
herence  squared  (y  2);  see  Eq.  (20).  Direct  comparisons  of  the  two  statistics,  y  5  and  p,  is  of  qual¬ 
itative  value  only  because  y2  is  frequency  dependent  while  px)  is  not. 
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The  seasonal  components  of  <5Z,6>Z)L.s;s-0  and  3LOD.m  (Fig.  4  on  page  23) appear 
to  be  two  in-phase  sinusoids  of  annual  and  semiannual  frequency.  The  linear  correlation. 
0.990,  is  significantly  higher  than  the  subseasonal  linear  correlation  of  0.853  (refer  to 
Table  1).  The  6LODvsso  has  slightly  greater  amplitude,  shown  by  its  greater  variance 
(Table  1).  Since  a  fundamental  assumption  (made  earlier)  is  that  the  d LOZ\sno  signal 
represents  the  total  of  all  atmospheric  (as  well  as  other,  unmodelled,  sources  of)  angular 
momentum,  it  is  not  surprising  that  the  seasonal  terms  of  the  atmosphere  are  less  than 
the  time  signals.  The  neglected  pressure  term  in  (13),  oceanic  currents,  and  stratospheric 
winds  are  probably  the  missing  inputs  to  balance  the  seasonal  terms.  Of  particular  note 
is  the  difference  in  amplitudes  found  in  mid  July  of  1984  and  1986.  The  Quasi- 
Biennual-Oscillation  (QBO)  peaked  at  exactly  those  two  times  (R.A.  Madden,  personal 
communication,  1988).  The  peak  was  for  easterly  stratospheric  winds;  enough  to 
qualitatively  account  for  the  significant  amplitude  difference  of  the  troposphere  and 
measured  LOD. 

Figs.  5  and  6  are  the  decompositions  of  the  seasonal  signals  depicted  in  Fig.  4 
on  page  23  into  annual  and  semiannual  parts.  Annual  terms,  Fig.  5  on  page  24,  appear 
smaller  than  those  determined  by  harmonic  analysis  (Table  3)  because  the  filter  transfer 
function  (required  to  separate  the  closely  spaced  annual  and  semiannual  terms)  atten¬ 
uated  the  annual  signal  65%.  Fig.  19  in  Appendix  C  is  the  plot  of  the  spectral  response 
of  the  filter  used  to  isolate  these  terms.  Residuals  of  dLOD^so  and  SLODUm  show  a 
small  biannual  bias  which  is  probably  associated  with  the  QBO. 

Fig.  6  on  page  24  contains  the  semiannual  part  of  the  seasonal  signal.  The 
spectral  response  of  the  filtering  (Fig.  20,  Appendix  C)  show  minimal  amplitude  losses 
due  to  the  filtering  process.  The  amplitude  from  the  plot,  about  0.3  ms,  agrees  well  with 
the  harmonic  analysis  (Table  3)  and  with  the  findings  of  Rosen  and  Salstcin  (1985). 
These  data  show,  as  other  data  have  in  previous  studies,  that  variations  in  the  length  of 
the  day  on  seasonal  timescales  (annual  and  semiannual  periods)  are  almost  entirely  due 
to  meteorological  influences. 

2.  Subseasonal  Components 

Subseasonal  is  defined  to  include  periods  from  20  to  150  days.  Fig.  7  on  page 
25  shows  the  close  match  of  the  subseasonal  signals  in  both  phase  and  amplitude.  The 
correlation  coefficient  of  the  two  series  is  0.853.  The  high  peaks  are  associated  with  pe¬ 
riods  around  50  days  while  the  lower  peaks  correspond  to  shorter  periods,  as  is  expected 
with  the/-1  power  spectrum  discussed  earlier.  The  50-day  peaks  are  perhaps  associated 


with  the  40-50  day  oscillations  in  the  tropics  first  described  by  Madden  and  Julian 
(1971).  The  average  amplitude  of  0.25  ms  is  slightly  larger  than  that  found  by  Eubanks 
et  al.  (1985)  for  the  period  1977-1981,  and  agree  well  with  the  findings  of  Morgan  et  al. 
(1985)  for  the  period  1981-1983. 

B.  GLOBAL  COMPARISONS  IN  THE  FREQUENCY  DOMAIN 

1.  The  Energy  Density  Spectra  of  SLODum  and  <5LODlsno 

Plots  of  SLOD^so  and  SLODtm  energy  density,  coherence  and  phase  are  con¬ 
tained  in  Fig.  12  on  page  28.  The  most  striking  feature  of  the  energy  density  plot  is  the 
large  peaks  at  annual  and  semiannual  periods.  This  is  consistent  with  the  largest  ampli¬ 
tude  oscillations  in  the  time  domain  (Fig.  2  on  page  22)  occurring  seasonally.  Also 
noteworthy  are  the  sharp  spikes  at  periods  of  13  and  seven  days  in  the  6L0DviSO  spec¬ 
trum.  These  are  high  frequency  tidal  oscillations  which  were  not  filtered  out  of  the  ori¬ 
ginal  time  data. 

The  slope  of  both  spectra  is  -2  (indicated  by  the  sloping  double  line  on  the  en¬ 
ergy  density  plot).  Eubanks  et  al.  (1985),  state  that,  while  the  cause  of  the/  2  dependence 
in  the  power  spectrum  is  not  clear,  it  shows  up  in  a  wide  variety  of  geophysical  phe¬ 
nomena.  Near  the  Nyquist  frequency  of  .5  day-1,  the  SLODvsso  spectrum  falls  off  at  a 
more  rapid  rate  than  expected  by  the / 2  relationship.  This  is  an  indication  of  the  lower 
noise  level  of  the  3LOD^%0  data.  Eubanks  et  al.  (1985)  experiments  were  conducted 
prior  to  recent  advances  in  more  accurate  astronomic  LOD  measurements;  thus  their 
data  probably  have  a  somewhat  noisier  spectrum  at  higher  frequencies  (personal  com¬ 
munication  with  J.O.  Dickey,  co-author  with  Eubanks  et  al.,  1985). 

At  periods  of  16  days  and  greater,  the  two  spectra  are  in  agreement.  In  addition 
to  the  broad  peaks  at  annual  and  semiannual  periods,  there  are  broad  peaks  at  75  and 
50  days,  and  narrower  peaks  at  33  and  25  days. 

2.  Coherence  Spectrum 

The  coherence  spectrum  in  Fig.  12  shows  strong  coherence  (or  correlation)  at 
seasonal  as  well  as  periods  centered  around  75,  50,  33,  25  and  20  days.  As  would  be 
expected,  these  periods  show  corresponding  power  in  the  energy  density  plots.  The  95% 
confidence  line  in  the  coherence  plot  is  derived  from  a  95%  significance  level  in  rejecting 
a  null  hypothesis  of  no  coherence  between  the  two  series.  The  coherence  plot  agrees 
well  with  that  of  Eubanks  et  al.  (1985).  Of  note  is  the  rapid  drop  off  in  coherence  at 
periods  around  125  days  (also  found  by  Eubanks  et  al.  1985).  This  may  be  related  to  a 
general  lack  of  power  in  both  the  SLODam  and  5LODcsno  spectra,  however,  the  loss  of 


power  does  not  completely  explain  the  steep  drop  in  coherence.  It  is  possible  that 
whatever  (weak)  oscillations  there  are  in  the  length  of  day  at  periods  around  125  days, 
they  are  not  of  meteorological  origin.  Significant  coherence  at  periods  as  low  as  20  days 
indicates  that  the  atmosphere  can  generate  short  period,  high  frequency,  changes  in  the 
LOD  .  Given  accurate  enough  time  and  wind  data,  fluctuations  of  even  one  or  two  day 
time  scales  will  eventually  be  resolved  (Lambeck,  1980). 

3.  Phase  Spectrum 

The  phase  plot  shows  that,  except  for  a  small  spike  at  125  days,  changes  of 
c5L<9Z)ltmand  <5LOZ)lsno  occur  within  five  days  of  each  other.  Generally,  phase  is  closest 
to  zero  where  coherence  is  highest,  indicating  changes  in  angular  momentum  are  rapidly 
transmitted  between  Earth  and  atmosphere.  Measurements  of  atmospheric  angular 
momentum  in  this  study  cannot  readily  be  decomposed  into  specific  torques  responsible 
for  transfer  of  momentum.  Swinbank  (1985)  found  that  shorter  period  fluctuations 
(from  daily  to  subseasonal)  result  from  mountain  torques  (pressure  differences  across 
north  south  oriented  mountain  ranges)  and  seasonal  variations  result  from  friction 
torques  (frictional  drag  at  the  surface  boundary'  layer  of  the  Earth  and  atmosphere). 

4.  Harmonic  Analysis 

Amplitude  and  phase  information  of  the  first  two  harmonic  frequencies  of 
5LOD^so  and  <5  LOD,m  were  calculated  by  year.  The  results  (Table  3)  compare  favora¬ 
bly  with  previous  studies  (Morgan  et  al.,  1985,  and  Rosen  and  Salstein,  1985),  although 
the  annual  amplitudes  found  in  this  study  show  more  variability  than  those  found  in  the 
previously  cited  studies.  Of  note  is  the  two-month  variation  in  phases  (Jan-Feb  for  an¬ 
nual  terms  and  May-Jun  for  semiannual  terms). 

C.  LATITUDE  COMPARISONS  OF  LOD 

Comparisons  of  SLODvs NO  with  tropical  and  extratropical  partitions  of  5LOD,m are 
shown  in  Figs.  13  through  16.  As  noted  above,  the  reason  for  computing  6LODtm  over 
specific  latitude  bands  (equatorial  and  extratropical  ,  ^gions,  respectively)  is  to  try  to 
identify  the  particular  geographical  region  that  is  contributing  to  the  high  coherence 
between  SLOD,tm  and  <5L(9Z?USN0.  Energy  density  plots  show  significant  drops  in  power 
at  periods  between  one  year  and  30  days  for  the  extratropical,  3LOD,tm-,s  and  3LODllm,0, 
scries.  The  two  tropical  series,  ( 6LOD,m[S  and  6LOD,m20)  show  power  losses  only  at  pe¬ 
riods  longer  than  100  days.  Coherence  for  the  extratropical  series  dropped  below  0.5  for 
most  periods  between  30  and  100  days;  for  the  tropical  series,  coherence  at  these  periods 
dropped  only  slightly.  Coherence  falls  off  rapidly  in  both  tropical  and  extratropical  series 


around  periods  of  125  days.  The  extratropical  series  fall  to  near  zero  at  these  periods. 
This  may  be  related  to  the  local  minimum  of  power  at  periods  from  80  to  125  days  in 
all  atmospheric  and  L’SN'O  energy  density  plots.  The  phase  plots  show  that  at  the  an¬ 
nual  period,  changes  in  the  atmosphere  tend  to  lead  changes  in  the  length  of  day  in 
extratropical  latitudes  (Fig.  14  and  especially  16)  while  they  tend  to  lag  changes  in  the 
length  of  day  in  the  equatorial  regions.  This  simply  reflects  the  fact  that  the  annual  cycle 
in  the  equatorial  region  is  shifted  or  delayed  somewhat  in  time  relative  to  the  middle 
latitudes. 

At  subseasonal  periods,  it  is  very  difficult  to  interpret  much  about  the  geographical 
origin  of  events  with  periods  shorter  than  about  30  days.  The  generally  high  power  and 
coherence  near  20-30  days  in  Fig.  12  does  not  seem  to  be  of  only  equatorial  or  only 
midlatitude  origin.  However,  the  evidence  very  clearly  indicates  that  the  tropics  and 
equatorial  regions  are  responsible  for  most  of  the  fluctuations  of  5LOD  in  the  40-to 
100-day  range.  Thus,  the  high  power  and  coherence  at  50-60  day  and  at  70-100  day 
periods  in  Fig.  12  also  appear  in  Figs.  13  and  15  (using  only  tropical  wind  data),  but  they 
are  completely  absent  in  Figs.  14  and  16  where  only  extratropical  wind  data  were  used 
in  computing  6LOD,lm.  The  dividing  line  between  the  tropical  and  extratropical  regions 
is  not  very  clear  because  only  two  latitude  bands  have  been  considered.  Also,  seasonal 
signals  are  still  apparent  in  the  tropical  regions,  and,  to  a  lesser  degree,  subseasonal 
(40-to  100-day  periods)  oscillations  are  seen  outside  the  tropics. 
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IV.  SUMMARY 

It  is  clear  that  changes  in  atmospheric  angular  momentum  can  account  for  the  vast 
majority  of  short  period  (one  year  or  less)  changes  in  the  LOD  .  Time  and  frequency 
domain  analyses  show  that  over  90%  of  the  variance  in  seasonal  and  subseasonal  peri¬ 
ods  of  5LODam  is  explained  by  5LODvsso.  Meteorological  activity  is  responsible  for  all 
changes  in  LOD  within  present  measurement  capabilities.  While  error  analysis  of  the 
atmospheric  and  time  series  is  beyond  the  scope  of  this  paper,  Morgan  et  al.  (1985)  ac¬ 
counted  for  most  of  the  variances  between  the  series  as  being  within  measurement  limits. 
In  addition,  Rosen  and  Salstein  (1985)  reduced  the  variance  due  to  unknown  sources  to 
less  than  0.04  (ms)2  for  annual  and  0.03  (ms)2  for  semiannual  terms  (measured  in  units 
of  amplitude).  Phase  differences  between  3LOD,lm  and  SLODlsso  are  very  near  zero,  in¬ 
dicating  that  an  almost  instantaneous  transfer  of  angular  momentum  occurs  between 
Earth  and  the  atmosphere.  Additionally,  when  the  pressure  term  is  included,  an  rms  er¬ 
ror  of  less  than  0.25  (ms)2  at  seasonal  periods  results  (Barnes  et  al.,  1983).  Other  studies 
(Eubanks  et  al.,  1985  and  1986),  conclude  that  measurement  errors  exceed  the  unex¬ 
plained  variance  of  6 LOD. 

The  value  of  establishing  the  strong  coherence  and  correlation  of  the  independently 
measured  signals  of  6LODlun  and  5L0DVSSCI  include: 

•  Astronomic  measures  of  LOD  can  be  used  as  an  independent  measure  of  the  state 
of  the  zonal  atmospheric  winds.  Analysis  and  forecasting  models  can  use  this  as  a 
check  of  data  consistency. 

•  Rapid  changes  of  LOD,  such  as  were  observed  in  the  El  Nino  event  of  1983,  can 
perhaps  be  used  to  infer  future  El  Nino  events  (Eubanks  et  al.  1986,  and  Rosen 
et  al.  1984).  Other  weather  phenomena  with  a  strong  zonal  component,  such  as  the 
40-50  day  tropical  oscillations,  can  be  temporally  located  by  real  time  filtering  and 
analysis  of  SLODlUn  data. 

•  Dickey  et  al.  (1986)  reported  that  atmospheric  angular  momentum  can  be  used  to 
infer  real  time  5  LOD.  This  is  required  for  extremely  precise  Earth  position  fixes  to 
navigate  spacecraft.  Real  time  calculations  of  atmospheric  angular  momentum  and 
5LOD,lrn  are  currently  computed  and  archived  by  NMC. 

The  findings  of  this  paper  are  in  agreement  with  previous  studies.  The  finer  resol¬ 
ution  of  the  L'SNO  time  data  beginning  in  1983  (the  BIH  time  data  also  had  similar 
improvements  starting  in  1983)  shows  that  atmospherically  and  astronomically  derived 
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3L0D  are  coherent  down  to  periods  of  20  days  or  less.  The  NOGAPS  wind  analyses 
compare  well  with  similar  studies  done  with  NMC  or  ECMWF. 

Comparing  latitude  bands  ofintegrated  5  LODtlm  with  5  LOD^%so  confirms  that  there 
is  spatial  variability  in  both  phase  and  amplitude  of  atmospheric  angular  momentum. 
Seasonal  changes  at  annual  periods  of  <5  LOD  appear  to  originate  primarily  in  the 
extratropical  regions.  Subseasonal  fluctuations  (30-100  day  periods)  in  3LOD,m have 
significantly  more  amplitude  (power)  and  coherence  with  <3LGOLsno  when  computed 
from  tropical  data  than  extratropical  data.  This  strongly  suggests  that  changes  in  the 
length  of  day  at  these  time  scales  is  produced  by  meteorological  events  in  the  tropics- 
e.g.  the  30-60  day  oscillations  of  Madden  and  Julian  (1971,  1972).  Further  study,  in¬ 
cluding  more  latitude  band  divisions  of  the  globe  is  required  to  be  more  exact  on  the 
position  of  the  dividing  line  between  the  tropics  and  extratropics.  This  study  indicates 
that  the  latitude  band  responsible  for  most  of  the  subseasonal  fluctuations  of  the  length 
of  the  day  is  approximately  20 °iV  to  20°S. 


Table  1.  LINEAR  CORRELATION  COEFFICIENTS 
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Table  3.  PHASE  AND  AMPLITUDES  OF  LOD 
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Note  1.  Amplitude  is  in  ms. 

Note  2.  Phase  is  the  peak  value  of  that  term  in  days  starting  with  Jan  1. 
For  the  semiannual  phase,  the  date  is  the  day  of  the  first  peak 
in  the  year. 

Note  3.  From  Rosen  and  Salstein  (1985).  Time  data  is  from  BIH.  NMC 
weather  data  extends  to  1  mb. 

Note  4  From  Morgan  et  al.  (1985).  Time  data  is  a  composite  of 

LLR,  VLB1,  and  BIH  estimates.  NMC  data  extends  to  100  mb. 


Fig.  1.  Unfiltered  time  series  of  6LOD^so  (top)and  5LODum  (bottom):  Note  the 
linear,  decreasing  trend  in  the  6LODvs^0  plot.  The  mean  has  not  been  re¬ 
moved  from  the  USNO  data. 
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Fig.  2.  Detrended  5LODMm  and  SLODvstio  series:  Top  line  is  USNO  data.  Bottom 
line  (atmospheric  data)  has  been  offset  by  -0.2  ms  for  clarity 


Fig.  3.  Smoothed  SLODltm  and  5LODvsyo  series:  A  23-point  moving  average  fil¬ 


ter  was  used  to  eliminate  noise  and  high  frequency  tidal  oscillations.  Heavy 
line  is  the  <5L0Z?USNO  series. 
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A  75-point  weighted  moving  average  filter  was  applied  twice  to  each  series; 
148  points  were  lost  from  the  ends. 


Fig.  5.  Annual  terms  of  t>LODvsyo  and  6LOD„\  Heavy  line  is  <5L0Z\SNOBottom 
graph  is  residual  of  dLOD^so  -  <5 LODtxm  shifted  downward  by  0.2  ms.  The 
slight  oscillation  of  this  line  is  at  a  period  of  2  years  and  is  probably  related 
to  the  Quasi- Biennual-Oscillation.  Frequency  transfer  function  is  Fig.  19 
on  page  48  of  Appendix  C. 


Fig.  6. 


Semiannual  terms  of  5LOD^0  and  SLODUm:  Heavy  line  is  SLOD^o 
Bottom  graph  is  residual  of  dLGDUSN0  -  5LOD,m  shifted  down  by  0.5  ms. 
Frequency  response  of  filter  is  Fig.  20  on  page  48  of  Appendix  C. 


Fig.  7.  Subseasonal  Oscillations  of  SLODxn  and  SLODvs^0:  Residual  of  series 
depicted  in  Fig.  4  on  page  23  subtracted  from  Fig.  3  on  page  23.  Bottom 
line  is  the  difference  of  the  two  series  shifted  down  0.5  ms. 


Fig.  8.  Comparison  of  SLODvayp^ and  5LODml)i  Heavy  line  is  <5L0Z\SNO 
seasonal  signal;  light  line  is  5LODumls  derived  from  integrating  atmospheric 
angular  momentum  within  15 0  of  the  equator. 


Fig.  11.  Comparison  of  6LODvsso.KMI>n^  and  5LODum-0:  Heavy  line  is  USNO  sea¬ 
sonal  signal  (from  Fig.  6);  light  line  is  atmospheric  signal  derived  from 
midlatitude  components  of  atmospheric  LOD  above  20°  N  and  below 
20°  S. 


Fig.  12.  Energy  Density,  Coherence,  and  Phase  Spectra  of  5LODMr,  and 
5LOD^iria:  Seven  point  spectral  smoothing  used  for  all  plots.  Solid  line 
in  the  power  spectral  density  spectral  plot  is  SLODvwo  and  dashed  line  is 
5LOD,m.  Dashed  lines  in  phase  plot  are  5  day  advances  (top  dashed  line) 
and  lags  (bottom  dashed  line)  of  <5 LODvs„o  with  respect  to  <5 LODxm. 
Horizontal  line  in  coherence  plot  represents  95%  confidence  that  the  two 
series  are  correlated. 
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Fig.  13. 


Energy  Density,  Coherence,  and  Phase  Spectra  of  5LODls>0  and 
8LODmiii  As  in  Fig.  12  except  the  atmospheric  signal  is  5LODtm  com¬ 
puted  within  15°  of  the  equator. 
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Fig.  16.  Energy  Density,  Coherence,  and  Phase  Spectra  of  SLODv 


8LODmkt  As  in  Fig.  12  except  the  atmospheric  signal  is  computed 
above  20°  N  and  below  20°  S. 
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APPENDIX  A.  FORTRAN  CODE  FOR  DECODING  NEDN  TAPES  AND 

COMPUTING  LOD 

//BENEDICT  JOB  (4507 ,0818) JAN  87',CLASS=G 
//*MAIN  0RG=NPGVM1. 4507P , LINES=20 

//  EXEC  FORTVCLG , PARM.  FORT=' 0PT(2) ,LVL( 66) 1 .REGION.  GO=2048K 
//FORT. SYSIN  DD  * 

C**************Vf**********I4SlP0RTANT,v**************************** 


ORIGINAL  AUTHOR  JAMES  BOYI  ,  NAVAL  POSTGRADUATE  SCHOOL  1986 
ADDITIONS  MADE  IN  1987  BY  W.  L.  BENEDICT 

THIS  IS  THE  MAIN  PROGRAM  TO  READ  THE  NEDN  TAPES  AND  GET  ZONALLY 
AVERAGED  WINDS  OUT.  BE  SURE  TO  CHECK  THE  INPUT  PARAMETERS  AT  THE 
END  OF  THE  PROGRAM  FOR  THE  DATES.  THE  OUTPUT  FILE  WILL  GO  TO  YOUR 
READER. 

************************* IMPORT AN t****************************** 
PGM  TO  READ  DATA  FROM  FNOC  NOGAPS  TAPE  GLOBAL  GRID  NEDN  FORMAT 

»»»»>»»»»»  LONG  RECORD  FORMAT  <<<<<<<<<<<<<<<<<<<<<<< 


REAL*4  DATA( 105 12) , D( 73 , 144) ,C( 144 , 73) 

REAL*4  Z0NAVG( 1 1 , 84 ) , C0SSQ( 7 3 ) , AVERGE , FACTRV( 8 ) , VERT( 73), 
&  LODATM(400) ,  OMEGA,  ISHELL,AVG,G,PI,LODSID,R 

INTEGER  I3UF(5284) 

INTEGER  IDATA( 10512) 

INTEGER*2  JDATA( 21024) 

INTEGER  IN 1 ( 24) ,KFTC( 100) 

INTEGER  BCD( 24) ,FTC( 100) 

INTEGER  IDFELD( 3,12), KOUNT( 6,12,6) 

INTEGER  DELTAU, COUNT, QQ,Q 
INTEGER  CALNDR( 740) 

REAL  FCTR( 12) ,DBASE( 12) 

INTEGER  ITAU0K(4) 

C 

EQUIVALENCE  ( IDATA( 1 ) , JDATA( 1 ) ) 

EQUIVALENCE  (DATA( 1) ,D( 1 , 1) ) 

C 

NAMELIST/ INPUT/  NSKPRC , 1X1 , 1X2 , IDF1 , IDF2 
C 

DATA  FCTR/  12*0. 01/ 

DATA  DBASE/ 12*0. 0/ 

DATA  ITAU0K/00, 12,24,36/ 

DATA  NUMTAU/01/ 

DATA  DELTAU/ 12/ 

DATA  FACTRV  /2*7500 , 10000, 3*5000 , 2*2500/ 

DATA  OMEGA  /7. 29E-5/ 

DATA  ISHELL  /7.04E37/ 

DATA  PI  /3. 1415926/ 

DATA  G  /9. 81/ 

DATA  LODSID  /86000. 0/ 
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DATA  R  /6370000. 0/ 

C 

READ( 05, INPUT) 

WRITE( 06, INPUT) 

C 

READ( 05 , 1275)  ( ( IDFELD( I , J) , 1=1 , 3) , J-IDF1 , IDF2) 

1275  FORMAT( 1 1( 3A1 , IX) ) 

WRITEC 06 , 1276)  ( ( IDFELD( I , J) , 1=1 , 3) , J=IDF1 , IDF2) 

1276  FORMATC IX, 'FIELDS  SELECTED  ' , 1 1( 3A1 , IX) ) 

C 

COUNT=0 

QQ=0 

NA=21135 

NB=0 

NCNT=0 

NFLE=0 

C 

C  SET  UP  ARRAY  FOR  CALENDAR 
CALL  KLNDR(CALNDR) 

C  SET  UP  COSINE  ARRAY  FOR  INTEGRATION 
DO  2300  1=1,  37 

COSSQ(I)  =  (SIN(2. 5*(I-l)*(2*PI/360. 0))**2) 

COSSQC  74-1)  =  COSSQ(I) 

2300  CONTINUE 
C 

DO  10  1=1,6 
DO  10  J=1 , 8 
DO  10  K=1 ,6 
10  K0UNT(I,J,K)=0 
C 

IDDHH1=IX1  -  ((IX1/10000)*10000) 

IDD1=IDDHH1/ 100 
IHH1«IDDHH1  -  IDD1*100 
C 

CALL  BUFFER( 10,1, 21135, IBUF) 

C  SKIP  TO  DESIRED  DATA  E.  G.  DEC  74  NSKPRC  =  806 
CALL  SKPREC( 10, NSKPRC, &99,&999) 

C 

C  IDFELD  IS  THE  NEDN  FIELD  IDENTIFIER  E.  G.  A01  =  SLP 
C  F00=500  Z,  T21  =  250MB  V, 

C 

C  1X1  &  1X2  ARE  THE  DTG  OF  FIRST  AND  LAST  PERIODS  (FORM  YYMMDDHH  IN  18) 
C 

C  *****************  TOP  OF  DATA  INPUT  LOOP  ************** 

55  CONTINUE 
C 

C  LONG  RECORD  FORMAT  JUST  USE  ONE  READ 
CALL  GETRECC 10,&99 ,&9999) 

NCNT=NCNT  +  1 
C 

CALL  GBYTE S( IBUF( 1) , IN1( 1 ) , 00 , 06 , 00 , 24) 

CALL  XBCD(24 , IN1( 1) ,BCD( 1) ) 

CALL  GBYTESC IBUF( 6 ) , IFTC ,28,4,0,0) 

NFTC=IFTC*4 

CALL  GBYTESC IBUF( 16) ,KFTC( 1) , 24 ,6 ,0 ,NFTC) 

CALL  XBCD(NFTC,KFTC( 1) ,FTC( 1) ) 
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CALL  XINT(BCD,5 ,4, ITAU) 

************  FORECAST  HOUR  (  TAU  )  CHECK  ************************ 

DO  610  K=1 .NUMTAU 
NNTAU=K 

IF( ITAU. EQ.  ITAUOK(K) )  GOTO  620 
610  CONTINUE 
GOTO  55 
620  CONTINUE 

****  CHECK  YYMMDDHH  TO  GET  THE  DATA  WE  NEED  ***** 

CALL  XINT( BCD ,9,8, I YMDH) 

IF  (IYMDH.LT.  1X1)  GO  TO  55 
IF  (I YMDH. GT. 1X2)  GO  TO  77 
IDDHHX=IYMDH  -  ( ( I YMDH/ 10000)* 10000) 

I DDX=I DDHHX/ 100 
IHHX=IDDHHX  -  IDDX*10Q 

NNDAY=(IDDX-IDD1)  +  ( IABS( IHHX-IHH1 )/12)  +  1 

*******  DATA  TYPE  CHECK  ********************** 

DO  640  K=IDF1,IDF2 

NNFLD=K 

DO  644  KK=1 , 3 

644  IF(BCD(KK)  . NE.  IDFELD(KK.K) )  GOTO  640 
GOTO  645 

640  CONTINUE 
GOTO  55 

645  CONTINUE 

***  DECODE  SCALING  INFORMATION,  GRID  SIZE 
CALL  GBYTESC IBUF( 6 ) , I UNIT , 8 , 5 , 0 , 0 ) 

CALL  GBYTES( IBUF( 6) .ISCLE, 13,6,0,0) 

CALL  GBYTESC IBUF(6) ,KEYBIT,26,2,0,0) 

CALL  GBYTESC IBUFC 7) ,M, 24, 12,0,0) 

CALL  GBYTESC IBUFC 8), N, 4, 12, 0,0) 

NR=M*N 

***  DECODE  DATA  ITSELF  AFTER  FINDING  WHERE  IT  STARTS. 

NBITS  =  C  21  +  IFTC)*24 
IWDSRT  =  CNBITS/32)  +  1 
ISKIP  =  NBITS  -  ( (NBITS/32)*32) 

NBTSDP=16 

CALL  GBYTESC IBUFC IWDSRT) , IDATA .ISKIP , NBTSDP , 0 , NR) 

CALL  GBYTESC IBUFC 27 ) , IDATA , 8 , NBTSDP , 0 ,NR) 

CDC  HAS  ONE'S  COMPLEMENT,  IBM  HAS  TWO'S  COMPLEMENT  FOR  NEGATIVE'S. 

DO  30  I  =  1 ,NR 

30  IF  ( IDATA(I)  .LT.  0)  IDATA(I)  =  IDATA(I)  +  1 

SET  UP  TO  UNPACK  16  BIT  DATA  POINTS.  THUS  JDATA  WHICH  IS  *2  HAS  ONE 
DATA  POINT  PER  WORD.  DATA  IS  IN  THE  RIGHTMOST  16  BITS  OF  IDATA  WHICH 
IS  *4. 

L=0 

NR2=NR*2 

IF  (ISCLE  .LT.  32)  GO  TO  40 


ISL=  2**( ISCLE-32) 

DO  17  I  =  2.NR2.2 
L=L+1 

IDUM=JDATA(I) 

17  DATA(L)  =  FLOAT( IDUM)/FLOAT( ISL) 

GO  TO  45 

40  ISL=2**( 32 -I SOLE) 

DO  18  1=2 ,NR2 , 2 
L=L+1 

IDUM=JDATA(  I) 

18  DATA( L)=FL0AT( IDUM)*FLOAT( ISL) 

45  CONTINUE 

C 

C  DBASE  =  5574.0  FCTR  =  0.01  FOR  500MB  F00 

C  DBASE  =  762.0  FCTR  =0.01  FOR  925MB  R00 

C  DBASE  =  0000. 0  FCTR  =  1.  00  FOR  SLP  A01 

C  DBASE  ~  0000.0  FCTR  =+0.01  FOR  V  AND  0.01  FOR  U. 

DO  88  KZ=1 ,NR 

88  DATA(KZ)=  DATA(KZ)*FCTR(NNFLD)  +  DBASE(NNFLD) 

C 

C  ORIENT  I,J  INDICES  SUCH  THAT  I=E-W  J=N-S 
C  REVERSE  SENSE  OF  J  SO  THAT  IT  INCREASES  NWD. 

DO  25  1=1,73 
11=73  -  (1-1) 

DO  25  J=1 , 144 
25  C( J, II)=D( I , J) 

C 

IHR=I YMDH- ( IYMDH/ 100 ) * 100 
IDAY=( ( I YMDH- ( IYMDH/ 10000)*10000) -IHR)/100 
C  **************  ZONAL  AVERAGING  LOOP*******************************’" 
C 

2999  DO  2115  J  =  1,  73 

AVG  =  0. 0 

DO  2116  1=1,  144 
AVG  =  C(I,J)  +  AVG 
2116  CONTINUE 

ZONAVG( NNFLD , J)  =  AVG/ 144. 

2115  CONTINUE 

IF  (NNFLD  .EQ.  11)  GOTO  3000 
GOTO  55 

3000  QQ  =  QQ  +1 

IF  (IYMDH  .NE.  CALNDR(QQ))  GOTO  3500 
GOTO  3600 

3500  WRITE(6,3501)  QQ,  CALNDR(QQ) .LODATM(QQ) 

3501  FORMAT( IX, 13 , IX, 18 , 1X.F9. 6 , 2X, ' DATA  MISSING  FROM  FNOC  TAPES') 
GOTO  3000 

3600  DO  2130  1=1,  73 

VERT(I)  =  ( ZONAVG( 1,1) +Z0NAVG( 2 , I ) ) *F ACTRV ( 1 ) 

&  +  ( ZONAVG( 2,1) +ZONAVG( 3,1) )*FACTRV( 2 ) 

6c  +  (ZONAVG(  3 , 1)+Z0NAVG(4, 1 )  )*FACTRV( 3) 

6c  +  (  ZONAVG(  4 , 1  )+ZONAVG(  5,1)  )*FACTRV(  4) 

6c  +  ( ZONAVG(5 , 1 )+ZONAVG( 6, 1)  )*FACTRV( 5) 

6c  +  ( ZONAVG(  6 , 1)+ZONAVG(  7 , 1)  )*FACTRV( 6) 

6c  +  (  ZONAVG(  7 , 1  )+ZONAVG(8 , 1)  )*FACTRV(  7) 
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&  +  ( Z0NAVG( 8, I)+ZONAVG( 9 , I) )*FACTRV(8) 

&  +  (ZONAVG( 7 , I)+Z0NAVG(8 , I) )*FACTRV(9) 

2130  CONTINUE 

APPLY  COS  LAT  TERMS 
DO  2140  1=1,  73 

2140  VERT(I)  =  VERT( I )*C0SSQ( I ) 

HORIZONTALLY  INTEGRATE 
AVERGE=0.  0 
DO  2145  1=2,72 

AVERGE  =  VERT(I)  +  AVERGE 
2145  CONTINUE 

AVERGE  =  AVERGE*  PI/72. 0 

CONSTANTS  OUTSIDE  THE  INTEGRAL  ARE  NOW  BROUGHT  IN: 

AVERGE  =  AVERGE*(2*PI*(R**3))/G 
CALCULATE  THE  CHANGE  IN  LOD 

LODATM(QQ)  =  1000. *(AVERGE*L0DSID/(0MEGA*ISHELL) ) 

WRITE(6, 3001)  QQ,  IYMDH,  LODATM(QQ) 

3001  F0RMAT( IX, 13 , IX, 18 , IX, F9. 6) 

GOTO  55 

99  CONTINUE 

77  WRITE(6 , 151)  IYMDH 

*********************  OUTPUT  LOOP  ************************** 
WRITE  (6,1020) 

WRITE  (6,1010)  (LODATM(Q),Q=l,QQ) 

1010  FORMAT( 7F10. 6) 

1020  F0RMAT( IX, 'THESE  ARE  THE  FINAL  LOD  VALUES  IN  SECONDS') 

151  F0RMAT(1X, ’PGM  ENDED  OK  IYMDH  =  ’,110) 

DO  835  1=1 ,NNDAY 
WRITE( 06 , 1235)  1,1X1 

1235  FORMAT( IX , / / , IX , ' NDAY  =  ’,I3,2X,'IX1  =  ',18) 

835  CONTINUE 

CALL  FILEND( IDF1 , IDF2) 

STOP 

**************  parity  ERROR  TRAP  *** ********************* 

9999  WRITE(Q6,120) 

120  FORMAT( IX, '  PARITY  ERROR  GETREC  ') 

STOP 

C  ****************  SKPREC  FAILURE  ************************ 

999  WRITE(06, 122) 

122  F0RMAT( IX , '  SKPREC  FAILURE  ') 

STOP 

998  CONTINUE 

WRITE( 06 , 1298)  NFLE 

1298  F0RMAT( IX, '  PGM  ENDED.  NUMBER  OF  FILES  PROCESSED  =  ',14) 
C  CALL  FILEND( IDF1 , IDF2) 

STOP 
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END 

&&&&&&&&&&&&&&&&&  SUBROUTINE  XINT  &&&&&&&&&&&&&&&&&&& 

SUBROUTINE  XINT(BCD, ISRT, NUM, IOUT) 

THIS  ROUTINE  TAKES  NUM  EBDIC  CHARACTERS  FROM  BCD 

STARTING  AT  ISRT  AND  CONVERTS  IT  INTO  A  INTEGER  NUMBER  PUT  INTO  IOUT. 
INTEGER  BCD( 24) , DIG ITS ( 11) 

INTEGER  INTGRS( 11) ,INTS(8) 

DATA  DIGITS/ ' 0 ' 1  , '  2 ' , '3' ,'4' ,'5' ,’6* ,'7' ,'8* ,'9* '/ 

DATA  INTGRS /0, 1,2, 3, 4, 5, 6, 7, 8, 9,0/ 

DATA  I BLANK/'  '/ 

IEND=ISRT  +  NUM  -  1 
DO  100  I=ISRT, IEND 
II  =  I  -  ISRT  +  1 
DO  200  J=l,ll 

IF(BCD(I).  NE.  DIGITS(J))GO  TO  200 
INTS(II)=INTGRS(J) 

GO  TO  100 
200  CONTINUE 

WRITE( 6, 1000) II, ISRT, NUM, IEND, (BCD(K) ,K=ISRT, IEND) 

1000  F0RMAT(2H  ,' ERROR, CHECK  YOUR  SUBROUTINE  XINT1 ,4(2X,I4) ,2X,24A1) 
STOP 

INTS(II)  =  0 
100  CONTINUE 
IOUT=00 
IFCTR=1 

DO  300  1=1, NUM 
II  =  NUM  -  I  +  1 
IOUT  =  IOUT  +  INTS(II)*IFCTR 
IFCTR  =  IFCTR*10 
300  CONTINUE 
RETURN 
END 

&&&&&&&&&&&&&&&&&  SUBROUTINE  XBCD  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 


SUBROUTINE  XBCD(N, INPUT, BCDOUT) 

SUBROUTINE  XBCD  CONVERTS  N  INPUT  CHARACTERS  FROM  CDC  EXTERNAL  BCD 
TO  IBM  EBCDIC.  THE  INPUT  CHARACTERS  ARE  STORED  IN  ARRAY  INPUT 
AND  THE  OUTPUT  (CONVERTED)  CHARACTERS  ARE  STORED  IN  ARRAY  BCDOUT. 


INTEGER  HEX(47)/Z31 , Z32 , Z33 ,Z34 ,Z35 ,Z36 , Z37 , Z38, Z39 ,Z21 ,222 ,Z23 , 

*  Z24,Z25 ,Z26 ,Z27 ,Z28 ,Z29 ,Z12,Z13,Z14,Z15,Z16,Z17,Z18,Z19, ZOA, ZO 1 , 

*  Z02 , Z03 , Z04 , Z05 , Z06 , Z07 , Z08 , Z09 , Z30 , Z20 , Z10 , Zll , Z3B , Z2B , Z2C , Z1B , 

*  Z1C,Z0B,Z0C/ , INPUT(N) 

REAL  CHAR(47)/ ' A  * , ' B ' ,'C' ,'D' .'E1 ,'F' ,'G' .'H' , *  I ' ,'j' ,'K' ,'L' , 

*  '  M'  ,'N'  ,'0'  ,'P'  ,  Q'  ,'R'  ,fS'  ,  T'  ,'U'  ,'V'  ,  W'  ,  X’  ,'Y'  ,  Z'  ,  V  ,  l'  , 

t/f  let  l^t  I  t  tgt  »rtt  »c»  If  I  I  f/l  »  I  lAl  1*1  In! 

9  »  /  »  •  »  V  »  t  )  9 

*  '@7, bcdout(n) 

DO  100  1=1, N 
DO  200  J=1 ,47 

IF( INPUT( I).  EQ.  HEX( J) )  GO  TO  150 
200  CONTINUE 


J=39 

150  BCDOUT(I)=CHAR(J) 


jw  v\rv\rvuw^j 


I 


$ 

; 


13 


100  CONTINUE 
RETURN 
END 

C  *********************  SUBROUTINE  KLNDR  ****************************** 
SUBROUTINE  KLNDR(CALNDR) 

INTEGER  CALNDR( 740) ,0 

C  JAN  (31  DAYS) 

CALNDR(l)  =  87010100 
DO  5001  1=2,62,2 

CALNDR( I )=CALNDR( I - 1 )  +  12 
CALNDR(  I+1)=CALNDR( I- 1)  +  100 

5001  CONTINUE 

C  FEB  (28  DAYS) 

CALNDR(63)  =  87020100 
DO  5002  1=64,118,2 

CALNDR( I)=CALNDR( 1-1)  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5002  CONTINUE 

C  MAR  (31  DAYS) 

CALNDR( 119)  =  87030100 
DO  5003  1=120,180,2 

CALNDR( I )=CALNDR( I - 1 )  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5003  CONTINUE 

C  APR  (30  DAYS) 

CALNDR( 181)  =  87040100 
DO  5004  1=182,240,2 

CALNDR( I)=CALNDR( 1-1)  +  12 
CALNDR(I+1)=CALNDR(I-1)  +  100 

5004  CONTINUE 

C  MAY  (31  DAYS) 

CALNDR( 241)  =  87050100 
DO  5005  1=242,302,2 

CALNDR( I)=CALNDR(I-1)  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5005  CONTINUE 

C  JUN  (30  DAYS) 

CALNDR( 303)  =  87060100 
DO  5006  1=304,362,2 

C ALNDR( I ) =C ALNDR( I - 1 )  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5006  CONTINUE 

C  JUL  (31  DAYS) 

CALNDR(363)  =  87070100 
DO  5007  1=364,424,2 

CALNDR( I)=CALNDR( 1-1)  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5007  CONTINUE 

C  AUG  (31  DAYS) 

CALNDR( 425)  =  87080100 
DO  5008  1=426,486,2 

CALNDR( I ) =CALNDR( I - 1 )  +12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5008  CONTINUE 

C  SEP  (30  DAYS) 

CALNDR(487)  =  87090100 


DO  5009  1=488,546,2 

CALNDR(I)=CALNDR(I-1)  +  12 
CALNDR( 1+1 )=CALNDR( I - 1 )  +  100 

5009  CONTINUE 

C  OCT  (31  DAYS) 

CALNDR(547)  =  87100100 
DO  5010  1=548,608,2 

CALNDR(I)=CALNDR(I-1)  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5010  CONTINUE 

C  NOV  (30  DAYS) 

CALNDR(609)  =  87110100 
DO  5011  1=610,668,2 

CALNDR( I)=CALNDR( 1-1)  +  12 
CALNDR(I+1)=CALNDR(I-1)  +  100 

5011  CONTINUE 

C  DEC  (31  DAYS) 

CALNDR(669)  =  87120100 
DO  5012  1=670,730,2 

CALNDR( I )=CALNDR( I - 1 )  +  12 
CALNDR( I+1)=CALNDR( 1-1)  +  100 

5012  CONTINUE 
RETURN 
END 

C*************************  SUBROUTINE  FILEND 
SUBROUTINE  FILEND( IDSRT, IDEND) 

C  KK1=( IDSRT  +  1)*10 

C  KK2=( IDEND  +  1)*10 

C  DO  778  KK=KK1 ,KK2 , 10 

C  DO  778  K=2 ,02 

C  IFILE  =  KK  +  (K  -  1) 

C  778  ENDFILE  IFILE 

C  ENDFILE  21 

C  ENDFILE  22 

C  ENDFILE  24 

C  ENDFILE  26 

RETURN 
END 


/* 

//LBCED.  SYSLIB  DD 


// 

DD 

// 

DD 

// 

DD 

// 

DD 

// 

DD 

UNIT=3 35  0 , VOL=SER=MVS 003, D I S  P=SHR , D SN=C 1012.  TAPEUT , 

// 

LABEL=( , , , IN) 

//*  NOTE  ONLY  ODD  NUMBERED  TAPES  HAVE  THE  OOZ  DATA  AND  120HR  FCSTS. 
//GO. FT10F001  DD  UNIT=3400-5 , VOL=SER=W16733 , 

//  DISP=( OLD, PASS) ,LABEL=( 1 ,NL, , IN) , 

/ /  DCB=( RECFM=FB ,BLKSIZE=21135, DEN=4 ) 

//GO.  SYSIN  DD  * 

& INPUT  NSKPRC=00 ,  1X1=87070100,  1X2=87123112,  IDF1=01,  IDF2=9&END 
C20  D20  E20  F20  G20  H20  120  J20  K20 
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APPENDIX  B.  FORTRAN  CODE  FOR  SPECTRAL  ANALYSIS 

C*********************************************************************** 


ENERGY  DENSITY,  CROSS  SPECTRUM,  PHASE,  AND  COHERENCE 
INPUT  FILES  ARE  UNITS  33,  34,  and  35 
OUTPUT  FILES  ARE  UNITS  41  THROUGH  45 

***********  *******************************  ***************************** 


PROGRAM  CRSPEC 
FARAMETERC  N'=2048 , M=1 1 ) 

DECLARE  ARRAYS. 

REAL*4  DATAX(N) ,  GX(N),  GXSUM(N/2  +  1) ,ATM1( 1622) ,ATM2( 1622) 

REAL*4  DATAY(N) ,  GY(N) ,  GYSUM(N/2  +  1) ,L0D( 1622) ,LIMPHl(N/2+l) 
REAL*4  DAY,XSUM(N/2+l) ,YSUM(N/2+l) 

REAL*4  CXY(N) ,  QXY(N) ,  PHASE(N) ,  COHER(N) ,  FREQ(N) ,LIMPH2(N/2+l) 
INTEGER*4  IWK(N+1) ,FILENO 

COMPLEX  AX(N) ,  AY(N),  GXY(N) ,  GXYSUM(N/2  +  1),  XYSUM(N/2+l) 

CHARACTER*60  TITLE 1 .TITLE2 .TITLE3 ,TITLE4,TITLE5 .TITLE6 ,TITLE7 
DEFINE  CONSTANTS. 

M  AND  N  ARE  DEFINED  IN  PARAMETER  SECTION;  N  IS  NUMBER  OF  SAMPLES 
M  IS  THE  POWER  OF  TWO  TO  GET  N 
XN  =  FLOAT(N) 

NF  =  N/2  +  1 

(NF  IS  INDEX  VALUE  FOR  NYQUIST  FREQ. ) 

DT  -  1 

(DT  IS  TIME  BETWEEN  SAMPLED  VALUES) 

T  *  FLOAT(N) 

(T  IS  TOTAL  TIME  FOR  SAMPLE;  HERE  IT  IS  EQUAL  TO  THE  SAMPLE) 

DF  =  1.  /T 

(DF  IS  DELTA  F,  FUNDAMENTAL  FREQ.) 

NS  =  1 

(NS  IS  NO.  OF  SAMPLES  TO  ENSEMBLE  AVERAGE) 

PI  =  3. 141593 
WCORR  =  . 875 

(WINDOW  CORRECTION) 

INITIALIZE  SUMMING  ARRAYS.  RETAIN  INFORMATION  UP  TO  THE 
NYQUIST  FREQ.  ONLY. 

DO  100  I  =  1,  NF 
GXSUM(I)  =  0. 

GYSUM(I)  =  0. 

GXYSUM(I)  =  (0. ,0. ) 

100  CONTINUE 
C 

c ***********************  main  loop  of  program  ************************** 

C  COMPUTE  ENERGY  DENSITY  SPECTRA  AND  ENERGY  DENSITY  CROSS  SPECTRUM 
C  FOR  EACH  SAMPLE  AND  ENSEMBLE  AVERAGE. 

C 

FILENO=45 
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DO  500  K  a  1,  NS 

:  READ  DATA. 

I  THIS  IS  FOR  15N/S  DATA 

DO  144  1=  1,  1622 

I  READ  (34,209)  ATM1(I) 

BELOW  IS  ONLY  FOR  20N  TO  20S  DATA. 

READ  (33,209)  ATM1(I) 

ATM1( I)  =  16. / 7  2 . *ATM1( I) 

READ  (35,200)  L0D( I) , ATM1( I ) 

ATM1( I )=ATM2( I ) -ATM1( I) 

144  CONTINUE 

200  FORMAT  (2(1X,F9.6)) 

202  FORMAT  (24X,F8.5) 

209  F0RMAT( 10X,F9.  6) 

TAPER  BOTH  ENDS  OF  THE  DATA 

CALL  TAPER( 1622 ,L0D) 

CALL  TAPER( 1622 , ATM1) 

ADD  ZEROS 

DO  152  1=1, N 

IF  ((I  .LE.  213)  .OR.  (I  .GT.  1835))  THEN 
DATAX( I)=0. 

DATAY( I)=0. 

ELSE 

DATAX( I)=LOD( 1-213) 

DATAY(I)=ATM1( 1-213) 

ENDIF 
152  CONTINUE 

WRITE  (37,101)  (DATAY(I) ,I=1,N) 

101  F0RMAT( 8F8.  4) 

TRANSFER  THE  DATA  TO  COMPLEX  ARRAYS  AX  AND  AY  FOR  FFT2C  INPUT. 
NOTE  THAT  IN  USING  THE  2ND  FORMULA  IN  FFT2C  DOCUMENTATION,  YOU 
MUST  INITIALLY  TAKE  THE  COMPLEX  CONJUGATE  OF  THE  DATA.  SINCE 
THE  IMAGINARY  PART  IS  ZERO,  THIS  STEP  IS  NOT  NECESSARY. 

DO  300  I  =  1,  N 

AX( I)  =  CMPLX  (DATAX(I),  0.) 

AY( I)  =  CMPLX  (DATAY(I),  0.) 

300  CONTINUE 

DO  FFT. 

CALL  FFT2C  (AX,  M,  IWK) 

CALL  FFT2C  (AY,  M,  IWK) 

MAKE  ADJUSTMENT  TO  GET  2ND  FORMULA. 

DO  325  I  =  1,  N 
AX( I)  =  CONJG  (AX( I) )  /  XN 
AY( I)  =  CONJG  (AY(I))  /  XN 
325  CONTINUE 

COMPUTE  THE  ENERGY  DENSITY  SPECTRA  GX(F)  AND  GY(F)  AND  THE 
CROSS  SPECTRUM  GXY(F)  FROM  FFT2C  OUTPUT  AX  AND  AY. 
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DO  350  I  =  1,  N 

GX( I )  =  CABS(AX(I)**2)*(2.  /DF) 

GY( I)  =  CABS( AY( I )**2)*( 2.  /DF) 

GXY(I)  =  AY( I )*CONJG( AX( I ) )*( 2.  /DF) 

350  CONTINUE 
C 

C  SINCE  YOU  MUST  ENSEMBLE  AVERAGE  THE  RESULTS  OVER  THE  NO.  OF  SAMPLES, 
C  ADD  GX( F) ,  GY(F) ,  GXY(F)  INTO  SUMMING  ARRAYS.  NOTE  THAT  EACH 
C  FREQUENCY  IS  SUMMED  INDEPENDENTLY. 

DO  400  I  =  1,  NF 
GXSUM(I)  =  GXSUM(I)  +  GX(I) 

GYSUM(I)  =  GYSUM(I)  +  GYfl) 

GXYSUM(I)  =  GXYSUM(I)  +  GXY(I) 

400  CONTINUE 
C 

500  CONTINUE 

C*********************  end  OF  MAIN  LOOP  ******************************** 
C 

C  AFTER  SUMMING  IN  ALL  OF  THE  ESTIMATES,  DIVIDE  BY  THE  NO.  OF 
C  SAMPLES  NS  AND  BY  THE  WINDOW  CORRECTION  FACTOR  WCORR. 

C  COMPUTE  THE  CO-SPECTRUM  CXY,  THE  QUAD-SPECTRUM  QXY,  THE 

C  PHASE  SPECTRUM,  AND  THE  COHERENCE  SPECTRUM.  CONVERT  THE 
C  PHASE  SPECTRUM  TO  DEGREES  FOR  PLOTTING. 

XNS  =  FLOAT(NS) 

C  MULTIPLY  BY  F**2 

DO  600  1=2,  NF 

GXSUM(I)  =  (GXSUM(I)  /  (XNS  *  WC0RR))*((FL0AT(I-1)*DF)**2) 
GYSUM(I)  =  (GYSUM(I)  /  (XNS  *  WCORR) )*( (FLOAT( I-1)*DF)**2) 
GXYSUM(I)  =  (GXYSUM(I)  /  (XNS  *  WCORR) )*( (FLOAT( I- 1)*DF)**2) 

600  CONTINUE 

C  APPLY  A  5  POINT  SMOOTHING  CURVE) 

DO  610  1=3 ,NF-2 

XSUM( I )=( GXSUM( I- 1 )+GXSUM( I)+GXSUM( I+1)+GXSUM( I-2)+GXSUM( 1+2) )/5. 
YSUM( I)=(GYSUM( I-1)+GYSUM( I)+GYSUM( I+1)+GYSUM( I-2)+GYSUM( 1+2) )/5. 
XYSUM( I )=( GXYSUM( I - 1 )+GXYSUM( I )+GXYSUM( 1+1 )+GXYSUM( 1-2) 

&  +GXYSUM( 1+2 ) ) / 5 . 

610  CONTINUE 

XSUM( 2 )=( GXSUM( 2 )+GXSUM( 3 ) ) /2. 

YSUM( 2 ) =( GYSUM( 2 ) +GYSUM( 3 ) ) / 2. 

XYSUM( 2 )=( GXYSUM( 2 )+GXYSUM( 3 ) ) / 2. 

XSUM( 3)=(GXSUM( 2)+GXSUM( 3)+GXSUM(4) )/3. 

YSUM( 3 )=( GYSUM( 2 )+GYSUM( 3 ) +GYSUM( 4 ) ) / 3 . 

XYSUM( 3 ) =( GXYSUM( 2 ) +GXYSUM( 3 ) +GXYSUM( 4 ) ) / 3 . 

C  RETURN  VALUES  TO  ARRAYS 
DO  676  1=1, NF 
GXSUM( I )=XSUM( I ) 

GYSUM( I )=YSUM( I ) 

GXYSUM( I )=XYSUM( I ) 

676  CONTINUE 


2  CREATE  FREQUENCY  ARRAY  AS  MULTIPLES  OF  THE  FUNDAMENTAL  FREQ.  DF. 

2  NOTE  THAT  THE  FIRST  BIN  IS  THE  MEAN  AND  THE  FREQ.  IS  ZERO. 

DO  700  I  =  1,  NF 
FREQ(I)  =  FLOAT(I-l)  *  DF 
700  CONTINUE 

DO  615  1=2, NF 

CXY(I)  =  REALCGXYSUM(I)) 

QXY(I)  =  AIMAG( GXYSUM( I ) ) 

IF  (CXY(I)  .  EQ.  0.)  THEN 
PHASE(I)  =  0. 

GOTO  679 
END  IF 

PHASE(I)  =  ATAN2(QXY( I) ,CXY( I))*( 180. /PI) 

IF  (PHASE(I)  . GE.  90.)  PHASE( I )=90. 

IF  (PKASE(I)  . LE. -90.  )  PHASEC I )  =  -90. 

679  IF  ( ( GXSUM( I )*GYSUM( I ) )  . EQ.  0.)  THEN 
COHER( I )=0. 

GOTO  615 
ENDIF 

COHER(I)  =  ((CXY(I)**2)+(QXY(I)**2))/(GXSUM(I)*GYSUM(I)) 

615  CONTINUE 
2  DIVIDE  BY  F**2 

DO  611  1=2, NF 

GXSUM(  I )  =GXSUM(  I )  /  ( ( FLOA'I'(  I  - 1 )  *DF )  **2  ) 

GYSUM( I ) =GYSUM( I ) / ( ( FLOAT( I - 1 ) *DF ) **2 ) 

GXYSUMC I )  =GXYSUM(  I )  /  (  (  FLOATC I  - 1  )*DF),V,V2  ) 

611  CONTINUE 

'i 

:  SET  LIMITS  FOR  PHASE  PLOT 

DAY=5. 0 
DO  833  1=2, NF 

LIMPHK  I )=FREQ(  I  )*DAY*360 
IF  (LIMPHl(I).  GE.  90.)  LIMPH1( I)=90. 

LIMPH2C  I  )=  -LIMPHK  I) 

833  CONTINUE 

2  PRINT  OUT  FREQ.,  GXSUM(F) ,  GYSUM(F) ,  CXY(F),  QXY(F) ,  PHASE(F) , 

:  COHER(F). 

1005  FORMAT( 14X, 1 X1 ,9X,'Y' ,9X,’C0-' ,6X,'QUAD') 

1010  FORMAT( IX, 'FREQUENCY' ,2X, ' SPECTRUM' , 3X  ' SPECTRUM' ,2X, ' SPECTRUM’ , 
&2X, ’SPECTRUM’ ,4X, ’PHASE’ , 2X, ' COHERENCE 1 ) 

1015  FORMAT( IX, 72( ' - ' ) ) 

:  WRITE( 36, 1005) 

:  WRITE( 36 , 1010) 

3  WRITE( 36 ,1015) 

DO  800  I  =  2,  4 

WRITE(FILEN0,810)FREQ( I) ,GXSUM( I) ,GYSUM( I) ,PHASE( I) , 

&  COHERC I) .LIMPH1C I) ,LIMPH2( I) 

800  CONTINUE 

DO  802  I  ■  5,  NF-5 ,  5 
DO  844  K=1 ,5 

WRITE(FILENO ,810)FREQ(K+I ) ,GXSUM( 1+3) ,GYSUM( 1+3) , PHASE( 1+3) , 

&  COHERC 1+3) , LIMPHK K+I ) ,LIMPH2(K+I) 

844  CONTINUE 
802  CONTINUE 
810  FORMATC  7F10.  5) 
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C  COMPUTE  THE  VARIANCE  BY  INTEGRATING  GXSUM(F)  AND  GYStJM(F) 

C  NUMERICALLY  BY  SUMMING  GXSUM(F)*DF  AND  GYSUM(F)*DF  OVER  THE 
C  APPROPRIATE  RANGE  OF  FREQUENCIES.  WRITE  OUT  THE  VARIANCES. 

SUM  =  0. 

SUMM  =  0. 

DO  900  I  =  1,  NF 
PRODX  =  GXSUM( I)*DF 
PRODY  =  GYSUM( I )*DF 
SUM  =  SUM  +  PRODX 
SUMM  =  SUMM  +  PRODY 
900  CONTINUE 

WRITE( FILENO ,910)  SUM 
WRITE( FILENO, 920)  SUMM 

910  F0RMAT(//,4X, 'VARIANCE  OF  LOD  SERIES  =  ' ,F9. 3) 

920  F0RMAT(//,4X,’ VARIANCE  OF  ATM  SERIES  =  ' ,F9. 3) 

C  COMPUTE  THE  CO-VARIANCE  BY  INTEGRATING  THE  CO-SPECTRUM  CXY(F). 
C  WRITE  OUT  THE  CO -VARIANCE. 

S  =  0. 

DO  2000  I  =  1,  NF 
PRODXY  =  CXY( I)*DF 
S  =  S+  PRODXY 
2000  CONTINUE 

WRITE (FILENO, 20 10)  S 
2010  F0RMAT(//,4X,' COVARIANCE  =  ' ,F9. 3) 

STOP 

END 

SUBROUTINE  TAPER(NTS ,X) 

C  MULTIPLYING  INPUT  DATA  BY  TAPER  FUNCTION 
DIMENSION  X(NTS) 

PI  =  3. 1415927 
XNTS  =  NTS 
MA  =  0. 1*XNTS 
MB  =  XNTS-MA 
DO  204  I  =  l.NTS 
IF( I.  LE. MA)  GO  TO  201 
IF( I. GE. MB)  GO  TO  202 
GO  TO  204 

202  XI=NTS-I 
GO  TO  203 

201  XI=I-1 

203  ARG  =  5. 0*PI*XI/XNTS 

X( I)  =  X(I)*(1. -COS(ARG)*COS(ARG)) 

204  CONTINUE 
RETURN 
END 

SUBROUTINE  M INMAX  (A, N, AMIN, AMAX) 

C 

DIMENSION  A( 1) 

C 

AMIN  =  A( 1) 

AMAX  =  A( 1) 

C 

DO  100  I  =  1,  N 
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APPENDIX  C.  FREQUENCY  RESPONSE  OF  DIGITAL  FILTERS 


Fig.  17.  Spectral  response  of  23-point  sin*  weighted  filter:  Lowpass  filter  designed 
primarily  to  remove  solid  body  tides  at  7  and  13.5  days  from  5£6>DUSNO 
data. 


Fig.  18.  Spectral  response  of  a  twice  run  75-point  sin*  filter:  Lowpass  design  per¬ 
mitted  removal  of  most  subseasonai  components.  Note  the  abscence  of 
side  lobes,  as  well  as  the  relative  steepness  achieved  by  running  the  filter 
twice. 


a 
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Fig.  19. 


Response  of  183-point  unweighted  moving  average  filter:  Used  to  isolate 
annual  terms.  A  weighted  filter  would  have  required  too  many-points,  and 
would  not  have  sufficient  vertical  steepness  to  include  annual  frequencies. 
Still,  almost  35%  of  the  annual  signal  is  attenuated  in  order  to  filter  out 
the  semiannual  terms. 


Fig.  20.  Filter  to  isolate  semiannual  terms:  Combination  of  twice  run  75-point 
sin5  weighted  filter  and  183-point  unweighted  subtraction  moving  average. 
Maximum  response  is  at  about  200  days. 
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